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Abstract. We use molecular dynamics computer simulations to calculate the fre- 
quency dependence of the specific heat of a Si02 melt. The ions interact with the 
BKS potential and the simulations are done in the NVE ensemble. We find that the 
frequency dependence of the specific heat shows qualitatively the same behavior as 
the one of structural quantities, in that at high frequencies a microscopic peak is ob- 
served and at low frequencies an a-peak, the location of which quickly moves to lower 
frequencies when the temperature is decreased. 

INTRODUCTION 

The dynamics of supercooled liquids can be studied by many different techniques, 
such as light and neutron scattering, dielectric measurements, NMR, or frequency 
dependent specific heat measurements, to name a few [1]. In order to arrive at 
a better understanding of these systems also various types of computer simula- 
tions have been used to supplement the experimental data. However, essentially 
all of these simulations have focussed on the investigation of static properties or 
have studied the time dependence of structural quantities, like the mean squared 
displacement of a tagged particle or the decay of the intermediate scattering func- 
tion. What these simulations have not addressed so far, apart from a noticeable 
study of Grest and Nagel [2] , is the time dependence of thermodynamic quantities, 
like the specific heat. The reason for the lack of simulations in this direction is 
that the accurate determination of this quantity in a simulation is very demand- 
ing in computer resources because of its collective nature. This fact is of course 
very regrettable since one of the simplest ways to determine the glass transition 
temperature in a real experiment is to measure the (static) specific heat. Using 
ac techniques it is today also possible to measure the frequency dependent specific 
heat, c(v), and thus to gain more insight into this observable [3]. What is so far not 
possible in real experiments is to measure c(v) at frequencies higher than 1MHz, 
and thus the influence of the microscopic dynamics, which is in the THz range, 



cannot be investigated. For computer simulations it is, however, no problem to 
study c(u) also at these high frequencies and in this paper we report the outcome 
of such an investigation for the strong glass former silica. 



MODEL AND DETAILS OF THE SIMULATION 

The silica model we use is the one proposed by van Beest et al. [4]. In this model 
the interaction <f>(rij) between two particles i and j a distance r\j apart is given by 
a two body potential of the form 

<P(r tJ ) = q -f^ + Aij eM-Bi^j) - ^ . (1) 

' ij ' ij 

The values of the partial charges qi and the constants A^, Bij can be found in 
Ref. [4]. Since the quantity we want to investigate, c(v), is a collective one, it is 
necessary to average it over many independent realizations. Thus the system sizes 
we used are rather small, 336 ions, despite the fact that the dynamics of such a small 
system will show appreciable finite size effects [5]. However, exploratory runs with 
larger systems showed that these effects do not change the results substantially. 
The simulations were done at constant volume using a box size of 16.8 A, thus at 
a density of 2.36g/cm 3 , close to the experimental value of the density, which is at 
2.2g/cm 3 . The equations of motion have been integrated with the velocity form 
of the Verlet algorithm with a time step of 1.6 fs. The temperatures investigated 
were 6100 K, 4700 K, 4000 K, 3580 K, 3250 K, and 3000 K. At all temperatures 
the system was first equilibrated for a time which is significantly longer than the 
typical a-relaxation time of the system at this temperature. 



RESULTS 



In real experiments the frequency dependent specific heat is usually measured in 
the NPT ensemble. Although algorithms exist with which the static equilibrium 
properties of a system can be measured in a simulation in this ensemble, these 
algorithms introduce an artificial dynamics of the particles and are therefore not 
suited to investigate the dynamical properties of the system in this ensemble. Hence 
we calculated the frequency dependent specific heat in the microcanonical ensemble. 
Whereas in the NPT ensemble the specific heat is related to the fluctuations of 
the enthalpy, in the NEV ensemble it is related to the fluctuations of the kinetic 
energy [2,6]. It can be shown that in this ensemble the specific heat at frequency 
v is given by 

~ 2/3 - K(t = 0) - il-nv It exp(i2nut)K(t) ' ^ 

where K(t) is the autocorrelation function of the kinetic energy E kin and is defined 

as 
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K(t) = ^ T - (E kin (t) -£ kin )(£kin(0) -E kiQ ) . (3) 

^kin L J 

Here E k i Q is the mean kinetic energy and iV is the total number of ions. The 
derivation of Eq. (2) can be found in Ref. [6]. 

From Eq. (2) we see that the relevant quantity is the autocorrelation function 
K{t). The time dependence of this quantity normalized by its value at time t = 0, 
is shown in Fig. 1 for all the temperatures we investigated. 

From this figure we recognize that for high temperatures K(t)/K(0) decays very 
quickly to a value around 0.13 and then goes to zero like a stretched exponential. 
With decreasing temperature the function shows a plateau at intermediate times, 
the length of which increases quickly with decreasing temperature. Such a time 
and temperature dependence is very similar to the one found for the relaxation 
behavior of structural quantities, such as the intermediate scattering function [7]. 
Apart from these features the curves for the lowest temperatures show also a local 
minimum at around 0.02 ps, the depth of which increases with decreasing tempera- 
tures. The existence of this dip, as well as the observed high frequency oscillations, 
can be understood by realizing that within the harmonic approximation, which 
will be valid at even lower temperatures, the correlator is closely related to the 
autocorrelation function of the velocity, which is well known to show such a dip. 

Using Eq. (2) we calculated from K(t) the frequency dependent specific heat c(u). 
The real and imaginary part of this quantity are shown in Fig. 2 for all temperatures 
investigated. Let us first discuss c'(v): For very high frequencies we expect this 
function to go to the ideal gas value of 1.5, since the configurational degrees of 
freedom are not able to take up energy at such high frequencies. With decreasing 
v the function shows a fast increase in the frequency range which corresponds 
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FIGURE 2. Frequency dependence of the real and imaginary part (times 2.0) of the specific 
heat for all temperatures investigated. The upper and lower dashed lines are the values of c' (v) 
for the harmonic solid and the ideal gas, respectively. 



to the microscopic vibrations. For low temperatures this regime is followed by a 
plateau the height of which corresponds to the static specific heat of the system if 
no relaxation would take place, i.e. to the specific heat of the vibrational degrees 
of freedom. Since, however, on the time scale of the a-relaxation time the system 
relaxes, c\v) shows at the corresponding frequencies a further upward step. This 
feature is related to the fact that for very long times, or small frequencies, those 
configurational degrees of freedom which are not of vibrational type are relaxing 
and thus can take up energy. At even smaller frequencies the curves then show a 
plateau, the height of which is the static specific heat of the system. We see that 
with decreasing temperature this value is decreasing but is always significantly 
above the harmonic value given by the Dulong-Petit value of 3.0, since the relaxing 
configurational degrees of freedom give rise to an enhancement of the static specific 
heat. 

We also note that the height of the first step in c'{y) [coming from the low 
frequency side] is the configurational part of the specific heat. We see that at the 
lowest temperature this height is rather small, 0.7/cs per particle, in agreement 
with the experimental observation for strong glass-formers. 

All these features can also be seen well in the imaginary part of civ). At high 
frequencies we have a microscopic peak which corresponds to the vibrational degrees 
of freedom. The location and the height of this peak is essentially independent of 
temperature. At intermediate and low temperatures a second peak is seen at low 
frequencies, the so-called a-peak. Its position depends strongly on temperature in 
agreement with the observation that the a-relaxation time of structural quantities 
increases quickly with decreasing temperatures [7]. From the location of this peak 



we can read off v max , the frequency scale of the relaxation of the specific heat. As it 
will be shown elsewhere [6], the product of z/ max with r(T), the a-relaxation time of 
the intermediate scattering function, is essentially independent of the temperature, 
thus showing the intimate connection between the frequency dependent specific 
heat and the structural relaxation, in agreement with the prediction of Gotze and 
Latz [8]. 

To summarize we can say that we have presented the results of a large scale 
molecular dynamics computer simulation of a realistic model of viscous silica to 
investigate the frequency dependence of the specific heat. In the frequency regime 
which is accessible also to experiments our results are in qualitative agreement 
with the experimental data [3]. At higher frequencies we see the influence of the 
vibrational degrees of freedom on c(u). Since no experimental data for c(u) is 
available for silica we are not able to compare the results of the present simulations 
with reality. However, in a previous investigation we have shown that the present 
model gives very good quantitative agreement of the static specific heat with the 
one of real silica [9] and thus it is not unreasonable to assume that the results of 
the simulation on the dynamic quantity is reliable also. 
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